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Abstract 

Analysis  of  periodic  pulse  trains  based  on  time  of  arrival  is  considered,  with  perhaps  very 
many  missing  observations  and  contaminated  data.  A  period  estimator  is  developed  based 
on  a  modified  Euclidean  algorithm.  This  algorithm  is  a  computationally  simple,  robust 
method  for  estimating  the  greatest  common  divisor  of  a  noisy  contaminated  data  set.  The 
resulting  estimate,  while  not  maximum  likelihood,  is  used  as  initialization  in  a  three-step 
algorithm  that  achieves  the  Cramer-Rao  bound  for  moderate  noise  levels,  as  shown  by 
comparing  Monte  Carlo  results  with  the  Cramer-Rao  bounds.  An  extension  using  multiple 
independent  data  records  is  also  developed  that  overcomes  high  levels  of  contamination. 
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1  Introduction 


Analysis  of  pulse  trains  is  a  long  standing  problem,  with  applications  including  radar 
[7,8,10,17,18],  communications  [8,9],  neurology  [1,5,22],  and  astronomy  [21].  For  example,  in 
radar  it  may  be  desired  to  accurately  estimate  the  pulse  repetition  interval  (PRI)  [19,25].  A 
similar  problem  is  frequency  estimation  in  additive  white  Gaussian  noise  at  moderate  to  high 
signal-to-noise  ratio  using  phase  data  [12,23].  At  its  most  fundamental,  pulse  train  analysis  is 
based  on  time  of  arrival  information  only,  as  might  be  obtained  from  a  matched  filter  or  other 
detector.  We  assume  time  is  highly  resolved  and  ignore  any  time  quantization  error.  In  this 
work  we  are  primarily  concerned  with  a  single  periodic  pulse  train  with  (perhaps  very  many) 
missing  observations  that  may  be  contaminated  with  outliers.  Our  data  model  for  this  case, 
in  terms  of  the  arrival  times  tj,  is  given  by 

tj  =  <t>  +  kjT  +  r)j,  j  =  1, . . . ,  N,  (1) 

where  T  is  the  unknown  period,  (j)  ~  W[0,T)  is  a  uniformly  distributed  random  phase,  the  kj' s 
are  positive  non-repeating  integers,  and  r]j  is  zero-mean  additive  white  Gaussian  noise  with 
variance  such  that  3cr„  <  We  let  S  be  a  sample  set  of  N  arrival  times,  i.e., 

S  =  {ti}U  (2) 

Without  loss  of  generality  we  assume  the  elements  of  S  are  indexed  in  descending  order,  so 
tj  >  tj+ 1  for  j  =  1, . . . ,  N  —  1.  This  assumption  will  be  useful  later  and  does  not  represent  any 
practical  difficulty.  This  model  allows  for  missing  observations  through  the  distribution  of  the 
kj’s,  e.g.,  with  N  even,  kj  =  N ,  N  —  2,  N  —  4, . . . ,  2  implies  observations  at  kj  =  N  —  1,  N  —  3, . . . 
are  missing.  Outliers  are  not  explicitly  included  in  (1),  although  they  are  considered  later.  The 
assumption  of  Gaussian  noise  can  be  relaxed  to  include  iid  non-Gaussian  noise,  although  the 
results  are  no  longer  optimal.  Later  the  assumption  on  4>  is  relaxed  to  allow  for  any  positive 
constant,  and  given  a  realization  of  S  we  regard  it  as  an  unknown  constant.  Equation  (1)  is 
also  useful  for  modeling  sinusoidal  zero-crossing  times  at  moderate  to  high  signal-to-noise  ratio 
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[20].  Alternative  data  models,  commonly  used  in  neurology,  assume  “integrate  and  fire”  pulse 
generation,  e.g.,  see  Brillinger  [5]  and  references  therein.  Cumulative  error  terms  may  also  be 
included  in  some  cases  [10]. 

Given  the  sample  set  S  the  problem  is  to  recover  the  period  T  and  possibly  the  phase 
< f> .  The  minimum  variance  unbiased  estimates  for  this  linear  regression  problem  take  a  least- 
squares  form  (see  equations  (3),  (4),  and  (7)).  However,  this  requires  knowledge  of  the  kj's. 
We  therefore  propose  a  multi-step  procedure  that  proceeds  by  (i)  estimating  T  directly,  (ii) 
estimating  the  kj'1  s,  and  (iii)  refining  the  estimate  of  T  using  the  estimated  fc/s  in  the  least- 
squares  solution.  This  estimate  is  shown  to  perform  well,  achieving  the  Cramer- Rao  bound  in 
many  cases,  despite  many  missing  observations  and  contaminated  data. 

The  direct  estimate  of  T  (step  (i)  above)  is  obtained  using  a  modified  Euclidean  algorithm 
[6].  This  approach  is  motivated  by  the  fact  that,  in  the  noise-free  case,  the  value  of  T  is  very 
likely  to  be  the  greatest  common  divisor  (gcd)  of  a  set  of  samples  {tj  —  for  N  >  10.  The 

modified  Euclidean  algorithm  is  a  computationally  simple,  robust  method  for  estimating  the 
gcd  from  a  noisy  contaminated  data  set. 

The  paper  is  organized  as  follows.  The  next  section  gives  the  maximum  likelihood  solution 
and  Cramer-Rao  bounds  for  estimating  T  and  (f>.  Our  analysis  has  led  us  to  work  with  the 
data  set  {tj  —  tj+ 1}  -  J^1  so  as  to  avoid  estimating  <fi,  which  can  be  unreliable.  For  completeness, 
related  work  based  on  a  point  process  (zero-one  time  series)  viewpoint  is  briefly  reviewed  in 
Section  3.  Next,  we  describe  our  modified  Euclidean  algorithms  for  estimating  T,  including  a 
robust  version  for  use  with  contaminated  data.  In  section  5  the  three-step  refined  estimation 
procedure  mentioned  above  is  developed.  Its  performance  is  compared  to  Cramer-Rao  bounds 
via  Monte  Carlo  simulation,  revealing  the  very  good  performance  of  the  algorithm  with  many 
missing  observations  and  contaminated  data.  Also,  an  approach  using  multiple  independent 
data  records  is  developed  that  overcomes  very  high  levels  of  contamination.  Before  concluding, 
we  mention  in  Section  6  the  “common  oscillator”  problem  and  its  solution  using  our  approach. 
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2  Maximum-Likelihood  Estimation 


Given  the  sample  data  set  S  from  (2)  we  may  write 
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where  kj  6  J\f  (the  natural  numbers)  and  kj  >  kj+i .  In  compact  form  this  is 


(3) 


t  =  X/3  +  77, 


(4) 


where  /3  =  [</>,  T]T  and  the  definitions  of  t,  rj ,  and  X  follow  from  (3).  We  eliminate  <j>  by 
forming  the  differences  yj  =  tj  —  tj+i  =  (kj  —  kj+1)T  +  (r)j  —  rjj+i),  yielding 
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where  Sj  =  r]j  —  rjj+i.  As  in  (4),  we  may  write  (5)  compactly  in  an  obvious  notation  as 


(5) 


y  =  XdT  +  6. 


(6) 


Using  (6)  relaxes  the  assumption  on  <f> ,  allowing  it  to  be  any  positive  constant. 

Equations  (4)  and  (6)  are  Unear  regression  problems  whose  least-squares  solutions  yield 
the  minimum- variance  unbiased  estimate  when  the  noise  is  zero-mean  Gaussian,  e.g.,  see  Kay 
[13].  Denoting  the  noise  covariance  matrix  as  Rv  =  E[t]T]t],  the  solution  to  (4)  corresponds  to 
maximum  Ukehhood  (ML)  estimation  and  takes  the  form  of  a  least-squares  estimate 

f3  =  (XtR?X)-'XtR;H.  (7) 

For  white  noise  (7)  simphfies  further  with  Rv  =  cr*I.  Similarly,  the  solution  to  (6)  is 

f  =  (XjRi'X^XjRj'y,  (8) 
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where  R g  =  E[66T\.  We  have  assumed  white  noise  so  that 

2<t*,  k  =  0 

E[6Jn+k]  =  k  =  ±  1  (9) 

0,  else 

k  / 

and  Rg  =  (r^Rg  where  Rg  has  2’s  on  the  main  diagonal,  — l’s  on  the  first  upper  and  lower 
diagonals,  and  zeros  elsewhere.  In  general  Rg  is  full  rank  and  its  inverse  can  be  expressed 
element-wise  as  j )  —  ij/(N  +  1)),  and  is  therefore  easily  computed  [12]. 

Note  also  that  knowledge  of  is  not  necessary  for  calculating  (7)  and  (8)  due  to  cancellation. 

Although  optimal,  use  of  these  estimators  requires  knowledge  of  the  coefficient  matrices 
X  or  Xd-  Obviously,  this  is  not  a  problem  if  there  are  no  missing  observations,  for  then 
kj  =  N  +  1  —  j  for  j  —  1, 2, ...  TV.  However,  when  observations  are  arbitrarily  missing  and  so 
the  k^ s  are  not  known  in  general,  one  is  faced  with  more  unknowns  than  equations  in  (4)  and 
(6).  One  possible  solution  to  this  dilemma  is  to  refine  the  estimate  of  T  through  a  multi-step 
procedure.  First  use  an  alternative  method  for  estimating  T  that  does  not  require  knowledge 
of  the  kj' s,  use  the  estimate  of  T  to  estimate  the  kj' s,  and  then  use  the  estimated  kj's  in  (4) 
or  (6)  in  an  attempt  to  refine  the  initial  estimate  of  T.  Later  we  take  this  approach  using  a 
modified  Euclidean  algorithm  to  obtain  the  initial  estimate  of  T  without  knowledge  of  the  kj's. 

Before  ending  this  section  we  describe  the  Cramer- Rao  bound  (CRB)  on  the  ML  estimates 
produced  by  (7)  and  (8).  The  pdf  of  the  observations  t  in  (1)  is  multivariate  Gaussian  given 

by 

p(t\/3)  =  (2 ir)-K/V,-N  exp  j^j(t  -  X/3)T(t  -  X/J)j  ,  (10) 

yielding  the  log-likelihood  function 

lnp(t|/3)  =  —~(ln  2tt  -  In  <rj)  -  ^(t  -  X(3)T(t  -  X(3).  (11) 

The  Fisher  information  matrix  is  found  to  be 

J  =  — R[V)g(V(glnp(t|/3))T]  =  -2 ;XTX ,  (12) 
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so  that  the  CRB  on  /?;  is 


var {fii  -  /?;}  >  [J  1]ii,  (13) 

where  /?»  is  the  ith  element  of  f3  and  is  the  (i,i)  element  of  J-1.  The  inverse  of  XT X  is 


readily  found  to  be 

(XTX)~1  = 

so  that 


Xj=i 

-  £f=1  k0 

N 

(14) 


and 


var{T  —  T}  > 


var{<^>  —  (f>}  > 


Ncl 


NT.uk2i-&Uk>y' 

Ef=1  k) 


(15) 


(16) 


NZf^k]-(  Ef=I%)2' 

From  (15)  we  observe  that  the  CRB  is  reduced  for  smaller  <r2.  Also,  for  fixed  JV,  it  is  reduced 
when  the  spread  of  the  kj’s  increases. 

From  the  model  of  (6)  we  have 

p(yl/3)  =  (2a)- exp  j~(t  -  Xtff  Rj\t  -  X^)|  ,  (17) 


with  as  =  2 <7V.  The  Fisher  information  matrix  for  this  case  is 

1 


J  =  -2XTd  Rj'Xa, 


with  the  CRB  given  by 


var  {T-f}><T26{XjR?Xd)-\ 
which  can  be  evaluated  numerically  given  Xd. 

Remark 


(18) 


(19) 


1.  We  may  form  ML  estimates  of  T  from  (7)  or  (8),  with  corresponding  CRB’s  given  by 
(15)  and  (19).  However,  we  must  be  careful  when  comparing  the  CRB’s  because  the  two 
estimates  are  based  on  different  data. 
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3  Analysis  of  Point  Processes 

Before  describing  the  modified  Euclidean  algorithm  (in  section  4)  we  briefly  discuss  some 
alternatives  to  the  least-squares  ML  estimates.  The  general  approach  is  to  view  arrival  times 

tj  as  generating  a  zero-one  time  series  or  delta-train  given  by 

N 

x(t)  =  £*(<-  tj)-  (20) 

j=l 

Using  (1)  we  generate  a  corresponding  x(t)  that  is  a  periodic  pulse  train  consisting  of  jittered 
data  with  missing  observations.  More  complicated  models  include  outliers  and  other  pulse 
trains  interleaved. 

In  searching  x(t)  for  the  presence  of  a  periodic  delta-train  s(t)  given  by 

*(<)  =  £«(*-  m,  (21) 

j=i 

it  is  natural  and  intuitively  appealing  to  circularly  convolve  x(t )  and  s(t).  This  is  equivalent 
to  (circularly)  applying  a  matched  filter  to  s(t).  Note  that  the  matched  filter  is  not  an  op¬ 
timal  detector  for  s(t)  because  the  noise  in  x(t )  manifests  itself  in  pulse-position  rather  than 
amplitude.  Because  we  are  dealing  with  zero-one  time  series  the  convolution  is  equivalent  to 
counting  pulses  for  a  specific  hypothesized  phase  (j)  and  period  T.  Indeed,  perhaps  the  concep¬ 
tually  simplest  approach  to  detecting  s(t)  is  to  use  histograms  [25].  While  counting  schemes 
are  straightforward,  they  can  be  cumbersome  to  implement,  and  the  resulting  histograms  may 
be  difficult  to  interpret  [25].  In  cases  of  multiple  interleaved  pulse  trains  more  sophisticated 
counting  schemes  rely  on  reducing  the  search  space  over  T  in  clever  ways  [17,18]. 

Another  approach  is  to  perform  spectral  analysis  of  the  point-process  x(t).  This  can  be 
accomplished  using  a  periodogram,  given  by 

i  " 

P-U)  =  .  (22) 

j=  1 

e.g.,  see  Bartlett  [2]  and  Brillinger  [4].  A  search  for  peaks  of  Px(f )  produces  candidate  estimates 
of  T.  Peaks  of  Px(f )  have  been  shown  to  yield  approximate  ML  estimates  of  T  with  respect  to 
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the  point  process  x(t)  [8].  Note  that  this  does  not  imply  that  use  of  Px(f )  is  ML  for  the  data 
set  S.  Let  the  noise  rjj  have  the  pdf 


f(v)  = 


gOL  COs(27T77/T) 

2^70(q:)  ’ 


(23) 


with  parameter  a,  where  Io(oc)  is  the  zero-th  order  Bessel  function  evaluated  at  a.  This  pdf 
was  used  by  Van  Trees  for  analyzing  phase-lock  synchronization  schemes  [24,  ch.  3].  Note  that 
f{r})  tends  to  a  uniform  pdf  as  a  — >  0,  and  tends  to  a  Gaussian  pdf  as  a  — >  oo.  This  is  a 
fortuitous  choice  for  the  noise  pdf,  as  it  results  in  a  maximum  likelihood  estimate  of  T  given 
by  (see  the  Appendix;) 

2 

(24) 

We  use  the  notation  Tx  so  as  not  to  confuse  it  (which  is  ML  for  (20))  with  other  period 
estimates,  e.g.,  the  ML  estimate  given  by  (7)  for  equation  (3)  corresponding  to  the  data  of  (1). 
It  may  similarly  be  shown  that  the  ML  estimate  of  <f>  for  this  pdf  is  given  by 

(25> 

These  estimators  depend  directly  on  the  exponential  sum  (as  does  the  periodogram),  and 
do  not  require  knowledge  of  the  pdf  parameter  a.  Use  of  (24)  requires  a  fine  search  over 
/  =  1/T.  This  may  result  in  a  significant  computational  load  because  the  FFT  is  not  a  good 
computational  option,  versus  direct  computation  of  the  exponential  sum,  given  a  large  number 
of  zeros  in  x(t)  [6].  Nevertheless,  with  no  other  information,  the  use  of  Px(f )  is  well  founded, 
especially  for  multiple  interleaved  pulse  trains  without  too  many  missing  observations. 

A  procedure  related  to  the  periodogram  is  the  use  of  circular  statistics,  which  is  a  method 
for  synchronized  averaging.  This  approach  has  been  adopted  for  neuronal  pulse  train  analysis 
[1],  as  well  as  for  PRI  detection  [7].  For  p  >  0,  a  circular  statistic  is  given  by 


Tx  —  argmaxr 


N 

E 

j=i 


e2nt'r 


(26) 
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where  [*J  is  the  floor  function.  Thus  {■)  is  the  fractional  part,  and  so  £  [0, 1 ) .  If  p  is 

equal  to  T/k  for  positive  integer  k ,  then  a  histogram  of  will  exhibit  a  peak.  The  function 

is  a  mapping  of  tj  onto  a  circle  (or  folding  onto  an  interval).  It  is  interesting  to  note  that 

E"^  =  Ee  (27) 

3= 1  3= 1 

so  that  the  Fourier  transform  of  x(t)/p  gives  information  similar  to  a  histogram  of  tj/ p  mapped 
onto  the  unit  circle  in  the  complex  plane. 

4  Modified  Euclidean  Algorithms  for  Estimating  T 

In  this  section  we  describe  modifications  to  the  Euclidean  algorithm  for  estimating  T  in 
(1)  that  do  not  require  knowledge  of  the  kj's  [6].  Without  knowledge  of  the  kj's  the  resulting 
estimates  are  not  ML.  However,  as  described  in  the  next  section,  the  estimates  produced  can 
be  refined  to  yield  results  comparable  to  ML  estimation  for  moderate  noise  levels.  Motivation 
comes  from  the  fact  that,  in  the  noise-free  case  with  (f)  =  0,  T  is  very  likely  the  greatest  common 
divisor  of  the  set  S  given  by  (2). 

The  Euclidean  algorithm  is  the  method  for  finding  the  gcd  of  a  set  of  positive  integers, 
e.g.,  see  Hardy  and  Wright  [11],  Leveque  [15],  or  Knuth  [14].  The  algorithm  is  a  division 
process  for  the  set  of  integers  Z.  It  is  based  on  the  property  that,  given  two  positive  integers 
a  and  b,  a  >  6,  there  exist  two  unique  positive  integers  q  and  r  such  that 

a  =  q-  b  +  r,0<r<b.  (28) 

If  r  =  0,  we  say  that  b  divides  a.  This  property  of  the  set  of  integers,  combined  with  the  fact  that 
if  a,  b  £  Z  \  {0},  then  a  •  b  /  0  (Z  has  no  zero  divisors),  make  Z  a  unique  factorization  domain. 
Therefore,  in  Z  every  non-zero  element  may  be  written  uniquely  as  the  product  of  powers  of 
irreducible  integers,  or  primes.  The  Euclidean  algorithm  involves  repeated  application  of  (28) 
until  the  remainder  is  zero.  This  yields  gcd(o,6),  which  is  the  product  of  the  prime  factors 
that  divide  both  a  and  b. 


9 


The  Euclidean  algorithm  for  finding  gcd(a,  b)  can  be  simply  stated  as  follows,  e.g.,  see 
Knuth  [14,  p.319]. 

1)  Find  remainder  of  a/b. 

2)  If  remainder  equals  zero,  done,  and  gcd(a,  b)  —  b. 

3)  Set  a  b,  b  <—  remainder,  go  to  step  1. 

We  demonstrate  the  algorithm  showing  that  gcd(252,198)  =  18  =  2  •  32. 

252  =  1-198  +  54 
198  =  3-54  +  36 
54  =  1-36  +  18 
36  =  2-18. 

The  algorithm  can  be  extended  to  find  the  gcd  of  a  set  of  more  than  two  integers  by 
finding  the  gcd  pairwise  [14,  p.  324].  We  can  also  compute  the  gcd  of  a  set  whose  elements 
are  multiples  of  a  fixed  real  number,  as  is  the  case  of  interest  here  given  the  set  S.  It  can  be 
shown  that 

gcd^T, . . . ,  kNT )  =  T gcd(&i, . . . ,  kN)  (29) 

(see  Leveque  [15,  pg.  16]).  This  more  general  notion  of  gcd  is  consistent  with  Euclid’s  original 
algorithm  [14,  pp.  317-320]. 

The  standard  Euclidean  algorithm,  involving  repeated  division,  is  intolerant  to  the  pres¬ 
ence  of  noise  in  (1).  We  have  therefore  developed  a  modified  Euclidean  algorithm  (MEA)  to 
gain  robustness  to  noise.  The  MEA  exploits  the  fact  that  [6] 

gcd(fci, . . . ,  kN)  =  gcd((&i  -  k2),  ( k2  -  k3 ), . . . ,  (fcjv-i  -  kN),  kN) .  (30) 

It  uses  repeated  subtraction  rather  than  division.  The  basic  iterative  MEA  is  stated  as  follows. 
Recall  that  S  is  assumed  to  be  in  descending  order. 
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Modified  Euclidean  Algorithm 


1)  After  the  first  iteration,  append  zero  to  the  data  set. 

2)  Form  the  new  set  with  elements  tj  —  tj+i. 

3)  Sort  in  descending  order. 

4)  Eliminate  elements  in  [0,770]  from  end  of  the  set  (770  defined  in  text  below). 

5)  Algorithm  is  done  if  left  with  a  single  element.  Call  it  to,  and  declare  T  —  t0.  Else,  go  to 
step  1. 

The  MEA  works  by  repeatedly  exploiting  (29)  and  (30).  The  subtraction  of  adjacent 
pairs  in  the  first  iteration  eliminates  the  phase,  and  results  in  the  simpler  data  form 

S'  =  ,  with  =  KjT  +  ,  (31) 

where  Kj  =  kj  —  kj+ 1  (see  remark  1  at  the  end  of  this  section).  After  the  first  iteration,  we 
maintain  the  minimum  in  the  data  after  differencing,  just  as  in  (30),  by  appending  zero  in  step 
1  before  differencing.  After  step  2  differences  that  are  separated  by  (approximately)  T  will 
yield  small  (ideally  zero)  values.  These  small  values  are  eliminated  from  the  data  in  step  5 
with  proper  selection  of  threshold  r]0.  Choice  of  770  is  dictated  by  the  distribution  of  the  noise. 

The  MEA  will  yield  good  results  for  small  number  of  data  samples  N.  We  have  shown 
that,  given  N  ( N  >  2)  randomly  chosen  positive  integers  {ki, . . .  ,fc/v}, 

P{gcd(ku. . . ,  kN)  =  1}  =  [((N))-1  ,  (32) 

where  P{-}  denotes  probability  and  C(n)  is  the  Riemann  Zeta  function  (Theorem  3.1  in  [6]). 
We  then  get  as  a  corollary 

JHgcdCff,, . . . ,  =  1}  =  [<(JV  -  I)]-1  ,  (33) 

(Corollary  3.1  in  [6]).  We  also  show  that 

[l-21-w]"1<[C(JV)]-1<l,  (34) 
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(Proposition  3.1  in  [6]).  Thus,  [£(AQ]  1  converges  to  1  from  below  very  quickly  as  N  increases. 
For  example,  [£(10)]  1  =  0.9990  and  [£(16)]  1  =  1.0000.  Combining  these  results  with  (29) 
implies  that 

gcd^T, . . . ,  Kn.-lT)  — ►  T  (35) 

with  high  probability  for  relatively  small  N,  converging  to  T  quickly  as  N  — ►  oo.  In  the 
noise-free  case  this  theory  tells  us  that  the  algorithm  will  yield  T  with  high  probability  from 
only  N  =  10  data  samples.  In  practice  the  algorithm  yields  a  good  estimate  of  T  for  small 
(N  —  10)  to  moderate  (N  =  100)  sample  sizes  [6].  We  note  that  the  convergence  rate  of  the 
iterative  MEA  depends  on  the  spread  of  the  s,  in  that  large  spread  requires  more  iterations 
before  the  differences  approach  T. 

Other  modifications  of  the  MEA  were  developed  and  tested  in  [6] .  One  modification  is  a 
single  iteration  algorithm  that  is  robust  to  the  presence  of  outliers.  The  single  iteration  version 
of  the  MEA  is  possible  when  the  data  set  S  is  such  that  kj  —  kj+i  —  1  occurs  for  relatively 
many  values  of  j.  Assuming  this  is  the  case  then,  after  step  3  in  the  first  iteration,  there  will 
be  a  cluster  of  data  that  is  symmetrically  distributed  about  T  (because  ij  —  tj+i  «  T  when 
kj  —  kj+ 1  =  1).  Thus,  a  single  iteration  algorithm  proceeds  by  identifying  the  cluster  of  data 
grouped  around  T  (after  step  3  above),  and  then  averaging  over  this  cluster  to  estimate  T.  The 
data  cluster  may  be  isolated  using  a  gradient  operator  and  data-adaptive  thresholding.  This 
approach  has  the  additional  advantage  of  being  robust  to  the  presence  of  arbitrary  outliers  in 
S.  This  robust  single  iteration  algorithm  is  summarized  as  follows. 

Robust  Single  Iteration  Algorithm 

1)  Given  the  set  S,  form  the  new  set  with  elements  tj  —  tj+i. 

2)  Sort  the  new  set  in  descending  order. 

3)  Apply  gradient  operator  and  obtain  T,  with  y0  =  0.6 T  (parameters  defined  in  text  below). 

4)  Obtain  T  by  averaging  the  data  over  the  lowest  step,  isolating  the  step  based  on  the 
lowest  two  jumps  of  height  y0  with  step- width  greater  than  x0. 
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This  algorithm  requires  three  parameters:  a  gradient  threshold  g0,  the  step-height  thresh¬ 
old  yo,  and  the  step- width  threshold  x<>-  The  gradient  is  estimated  in  step  3,  with  large  gradient 
values  indicating  a  step  or  “edge”  in  the  data.  We  have  employed  a  simple  gradient  estimator 
by  convolving  with  an  impulse  response  given  by  [—1, 0, 1].  A  data-adaptive  gradient  threshold 
g0  is  selected  as  25%  of  the  maximum  gradient  value,  and  data  points  above  this  threshold  are 
assumed  to  correspond  to  the  step  edges.  The  choice  of  25%  is  based  on  experience  gained  by 
extensive  numerical  results,  and  may  typically  be  halved  or  doubled  without  greatly  affecting 
the  results.  The  rightmost  (i.e.,  smallest  in  average  amplitude)  cluster  in  the  data  between 
edges  corresponds  to  the  cluster  of  data  around  T.  We  therefore  take  the  rightmost  step, 
whose  height  we  denote  T,  as  a  coarse  estimate  of  T.  We  then  use  T  to  set  the  value  of  y0 
via  t/o  ~  0.6  T.  The  desired  data  cluster  is  then  obtained  as  the  smallest  (rightmost)  cluster 
between  two  edges  with  height  y0  such  that  the  number  of  data  points  in  the  cluster  is  greater 
than  xo- 

The  use  of  x0  represents  a  tradeoff  between  robustness  to  noise  and  outliers  versus  the 
possibility  that  the  algorithm  will  not  yield  a  solution.  For  higher  values  of  x0,  the  single 
iteration  algorithm  may  not  always  yield  a  solution,  but  will  generally  yield  reliable  estimates 
when  it  does.  For  lower  values  of  x0,  in  the  presence  of  many  outliers,  the  algorithm  will  likely 
always  provide  a  solution  but  may  yield  an  estimate  that  is  not  close  to  the  true  T.  In  our 
experience,  these  incorrect  estimates  are  often  less  than  the  true  T,  caused  by  a  falsely  identified 
data  cluster  made  up  of  outliers.  Thus  the  number  of  data  samples  within  the  selected  cluster 
is  a  good  measure  of  the  reliability  of  the  resulting  estimate  T. 

Remarks 

1.  It  is  easily  possible  to  generate  low  noise  or  even  noise- free  data  sets  such  that  the  MEA 
does  not  yield  the  true  period,  even  for  an  arbitrarily  large  number  of  data  points.  For 
example,  consider  ki  =  N,  k2  =  N  —  2,  ks  =  N  —  4, . . .  with  true  underlying  period  T  =  1. 
To  eliminate  the  phase,  we  do  not  save  the  first  minimum.  For  this  data,  the  MEA  will 
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give  T  —  2.  However,  for  even  10  data  points,  Theorem  3.1  and  Corollary  3.2  in  [6]  say 
that  such  data  sets  are  extremely  unlikely  to  occur.  Another  way  of  saying  this  is  that  if 
a  large  set  of  low  noise  or  noise-free  data  leads  the  MEA  to  an  incorrect  solution,  then 
the  corresponding  kj’s  exhibit  a  consistent  pattern. 

2.  The  computational  load  of  the  robust  single  iteration  algorithm  is  quite  small,  being 
roughly  2N  additions  plus  NlogN  operations  for  a  fast  sort.  This  assumes  convolution 
with  the  gradient  operator  [—1,0,1],  which  can  be  accomplished  with  N  additions  only. 

5  Achieving  ML  Performance  with  Missing  Observations  and  Out¬ 
liers 

The  MEA’s  described  in  the  previous  section  provide  low  complexity  methods  for  esti¬ 
mating  T  that  are  robust  and  perform  reasonably  well,  but  fall  short  of  the  CRB.  To  improve 
on  the  MEA  estimate  of  T  it  is  desired  to  estimate  the  coefficient  matrices  X  or  Xd  in  (4) 
or  (6),  respectively.  The  differenced-data  equation  (6)  is  preferred  over  equation  (4)  due  to 
the  difficulty  of  estimating  the  phase  <f>  (see  the  remark  below).  If  T  were  known  then  Xd 
could  be  estimated  using  (1  /T)  y.  Ideally,  this  estimate  is  composed  of  positive  integers,  but 
imperfect  knowledge  of  T  and  the  presence  of  noise  will  generally  yield  an  estimate  of  Xd  that 
has  non-integer  components.  We  therefore  propose  to  estimate  Xd  via 

(36) 

where  Tmea  is  the  estimate  of  T  obtained  via  the  MEA,  and  round]®]  =  \x  +  |J  denotes 
rounding  to  the  nearest  integer.  A  refined  estimate  of  T  is  then  obtained  by  using  Xd  in  (8) 
yielding 

T  =  {XjRJ,Xd)-1XjRJ,y-  (37) 

This  result  approaches  the  optimal  minimum  variance  performance  when  Xd  is  close  to  Xd- 
The  refinement  algorithm  is  summarized  as  follows. 


Xd  =  round 


[Tmea 


Refined  Estimation  Algorithm 


1)  Estimate  T  via  MEA  (  see  section  4),  calling  this  estimate  Tmea- 

2)  Estimate  Xj  via  (36). 

3)  Refine  the  estimate  of  T  using  Xj  in  (37),  calling  this  estimate  T. 

Performance  analysis  of  the  estimate  Tmea  depends  not  only  on  the  distribution  of  the 
noise  r/j,  but  also  on  the  distribution  of  the  kj’s.  Performance  analysis  of  the  iterative  MEA 
is  complicated  because  the  analysis  involves  order  statistics  and,  after  the  first  iteration,  the 
noise  is  no  longer  iid  [6].  We  focus  here  on  the  single  iteration  MEA,  and  utilize  an  iid  Bernoulli 
process  to  model  the  kj's. 

Suppose  tj  =  <f>  +  kjT  +  Tjj  is  an  element  of  S  for  some  kj.  Then,  the  probability  that 
kj+ 1  =  kj  —  1  is,  or  is  not,  an  element  of  S  is  modeled  by 

Prob(&j+1  =  kj  —  1)  =  1  —  A 

Prob^j+i  ^  kj  —  1)  =  A,  (38) 

where  0  <  A  <  1  is  the  Bernoulli  parameter.  If  kj  —  1  is  not  accepted  to  the  set,  then  kj  —  2 
is  accepted  with  probability  1  —  A  or  rejected  with  probability  A,  and  so  on.  Assuming  that 
ideally  kj  =  N  +  1  —  j,  j  —  1, 2, . . . ,  N  then  (38)  is  regarded  as  a  missing  observations  model. 
The  Bernoulli  model  is  commonly  employed  for  missing  observations  in  time  series  analysis, 
e.g.,  see  [3].  For  example,  with  A  =  0.25,  we  expect  25%  of  the  possible  observations  to  be 
missing.  It  follows  that,  after  a  single  iteration  of  the  MEA,  we  expect  (1  —  A )(N  —  1)  data 
samples  to  be  symmetrically  distributed  around  the  true  value  of  T.  Recall  that  these  samples 
arise  when  tj  —  tj+i  =  (kj  —  kj+i)T  +  (rjj  —  r]j+i)  is  such  that  kj  —  kj+i  =  1.  The  single  iteration 
MEA  estimates  T  by  identifying  and  averaging  the  data  cluster  around  T.  Under  the  white 
Gaussian  noise  assumption  we  therefore  have 

E[Tmea 3  =  T,  (39) 
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with  variance 

w{f“j}  =  (l-Apr-i)'  (40) 

Analysis  of  T  produced  by  (37)  is  difficult  due  to  the  nonlinearity  introduced  when 
estimating  Xd  in  (36).  We  therefore  quantify  the  behavior  of  T  by  comparing  Monte  Carlo 
simulation  with  theoretical  CR  bounds.  As  can  be  expected,  the  nonlinearity  of  (36)  leads  to  a 
threshold  effect  in  the  behavior  of  T  that  will  be  evident  when  comparing  the  sample  variance 
of  T  with  the  CRB. 

Data  without  outliers 

First  we  consider  noisy  data  without  outliers.  The  data  set  S  was  generated  with  N  =  100 
data  samples  as  in  (2),  where  the  fcj’s  were  generated  according  to  the  Bernoulli  model  of  (38) 
with  parameter  A.  The  phase  <f>  was  selected  randomly  from  [0,  T )  for  each  realization.  Without 
loss  of  generality  we  set  T  =  1  in  all  examples.  The  noise  rjj  was  white  Gaussian  with  variance 
cr^.  In  order  to  express  the  effect  of  the  noise  as  a  percentage  jitter  in  tj  we  use  the  definition 

%  jitter  =  )  x  100  (41) 

and  consider  av  values  such  that  0  <  %  jitter  <  50%.  In  all  cases  200  Monte  Carlo  trials  were 
conducted  for  each  value  of  %  jitter.  Results  are  plotted  for  a  normalized  mean-square  error, 
given  by 

(  T 2 

normalized  mean-square  error  =  10  log  -  — — 

H  6  \MSE 

where  MSE  was  obtained  experimentally  comparing  T  with  the  known  true  value  T.  The  MEA 

used  in  all  cases  was  the  robust  single  iteration  algorithm  described  in  Section  4. 

Figure  1  compares  three  estimates  of  T  with  the  CRB,  with  A  =  0.25  (25%  missing 
observations)  and  MEA  parameter  xq  —  5.  The  CRB  was  obtained  from  (19).  Note  that,  using 
the  normalization  of  (42),  the  CRB  lies  above  the  estimates  in  all  our  figures.  Performance 
of  Tmea  is  always  worse  than  the  CRB,  due  to  a  lack  of  knowledge  of  the  kj: s.  Two  refined 
estimates  are  also  shown.  Refinement  based  on  (36)  and  (37)  performs  well,  achieving  the 
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CRB  for  up  to  20%  jitter  and  thereafter  exhibiting  a  threshold  effect,  ultimately  degrading 
back  near  the  performance  of  Tmea • 

For  comparison,  performance  of  a  refinement  algorithm  based  on  (7)  is  also  shown  in 
figure  1.  We  estimate  X  using 

X  -=  round  — - (t  —  $x)  ,  (43) 

Tmea 

where  <f)x  is  the  phase  estimate  of  (25).  The  estimate  X  was  then  used  in  (7)  to  simultaneously 
refine  Tmea  and  (j>x .  This  approach  yields  relatively  poor  performance  due  to  the  poor  quality 
of  (f>x.  For  higher  values  of  A  we  have  observed  little  or  no  improvement  over  Tmea  for  even 
low  values  of  %  jitter,  so  further  experiments  are  not  reported  for  this  approach. 

Experimental  values  corresponding  to  figure  1  are  tabulated  in  table  1.  Experimental 
mean  and  standard  deviation  (std)  of  Tmea  and  the  refined  estimate  T  are  shown.  The  CRB 
is  given  in  terms  of  its  square  root  for  comparison  with  the  experimental  std’s.  The  std  of  T  is 
quite  close  to  the  CRB  for  up  to  20%  jitter,  while  stdl?^^}  is  one  to  two  orders  of  magnitude 
worse,  showing  the  level  of  improvement  achieved  by  refinement  of  Tmea •  Note  however  that, 
depending  on  the  application,  Tmea  is  a  reasonable  estimate  in  and  of  itself.  For  example,  at 
35%  jitter  std  {Tmea}  is  ~  1%  of  T  =  1.  Extensive  numerical  results  for  Tmea  are  given  in  [6]. 

Figures  2  and  3  show  results  for  A  =  0.5  and  A  =  0.75,  respectively.  All  other  parameters 
match  those  of  figure  1.  Beyond  A  =  0.5  there  is  a  lowering  of  the  threshold  for  which  T  achieves 
the  CRB,  as  might  be  expected  given  the  sparse  nature  of  the  data.  Note  the  change  in  CRB 
with  changing  A,  such  that  better  performance  is  possible  for  higher  A.  This  is  explained  by 
noting  that  for  a  fixed  number  of  data  samples  N,  a  higher  value  of  A  implies  a  greater  spread 
of  the  kj's  (see  the  discussion  following  the  CRB  equations  (15)  and  (16)). 

Data  with  outliers 

Further  Monte  Carlo  experiments  were  conducted  with  outliers  introduced  into  the  data. 
The  data  set  S  was  generated  as  before  and  then  arbitrary  outliers  were  added  by  selecting 
from  a  uniform  distribution  over  the  range  of  (outlier-free)  S,  i.e.,  over  [min(S),  max(5)].  For 
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example,  with  N  —  100  and  5%  outliers,  S  consists  of  95  data  samples  obeying  (1),  and  5 
arbitrary  samples. 

Naturally,  we  wish  to  somehow  remove  the  outliers  from  the  original  data.  However,  with 
missing  observations,  noise,  and  arbitrarily  placed  outliers,  an  outlier  removal  strategy  is  not 
obvious.  We  propose  using  the  MEA,  because  its  robustness  allows  it  to  function  in  this  case. 
So,  instead  of  reusing  all  the  data  in  step  3  of  the  refined  estimation  algorithm,  we  instead 
keep  only  the  data  that  is  utilized  in  forming  Tmea •  Let  the  subset  S'  C  S  denote  this  data, 
with  elements  i,-,  j  =  1, . . .  ,7V7  <  N.  The  reduced  set  S'  can  be  written  in  the  matrix-vector 
form  of  (3)  using  y-  =  t'-  —  i'-+1,  which  gives 


y'  =  x'dfMEA  +  s'. 


(44) 


Following  (36)  we  form 


Xj  =  round 


1 


(45) 


[Tmea 

and  apply  this  to  (44),  whose  solution  is  of  the  same  form  as  (37).  We  denote  this  refined 
estimate  based  on  S'  as  T'.  Note  that  Rj1  is  the  same  as  before,  except  for  a  size  change. 

Figure  4  compares  using  S  versus  S'  in  the  refinement  algorithm.  Here,  A  =  0.25,  Xo  —  15, 
and  the  data  was  contaminated  with  5%  outliers.  The  reuse  of  the  full  data  set  S  does  not  yield 
significant  improvement  over  Tmea •  Note  that  T  is  actually  worse  than  Tmea  for  jitter  <  8%. 
This  is  because  the  kj  values  corresponding  to  outliers  are  not  integers,  but  are  forced  to  be 
so  when  forming  Xj.  As  the  noise  is  reduced  this  induced  rounding  error  becomes  dominant. 
Also  shown  in  figure  4  is  the  normalized  MSE  of  T',  whose  performance  mimics  outlier-free 
results,  in  that  a  threshold  effect  is  evident  above  15%  jitter.  Below  the  threshold  T'  very 
nearly  achieves  the  CRB. 

That  T'  does  not  completely  achieve  the  CRB  is  due  to  the  inclusion  of  arbitrary  outliers 
in  S',  which  appears  to  be  generally  unavoidable.  Let  sm  £  S'  denote  the  rath  outlier  in  S, 
and  consider  two  cases.  First,  if  two  outliers  are  such  that  sm  —  sn  &  T,  then  they  will  be 
accepted  by  the  MEA  as  good  data.  Second,  if  sm  &  kT  +  (f>  for  some  integer  &,  then  sm  will  be 
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treated  as  valid  by  the  MEA.  The  likelihood  of  such  cases  obviously  increases  with  the  number 
of  outliers  present.  The  first  case  will  generally  have  more  of  a  negative  impact  on  T'  due  to 
contamination  of  Xd.  It  is  unlikely  that  outliers  in  either  case  can  be  identified  and  removed 
without  further  information.  However,  if  other  information  is  available,  such  as  knowledge 
of  their  distribution,  then  exploiting  this  to  remove  them  is  likely  to  be  both  possible  and 
worthwhile. 

An  examination  of  T'  reveals  interesting  behavior  near  the  threshold  (around  15%  jitter 
in  figure  4).  Plots  of  T'  are  shown  in  figures  5  and  6,  with  the  %  jitter  fixed  for  each  plot. 
Figure  5  shows  T'  over  200  Monte  Carlo  runs  for  15%  jitter,  exhibiting  low  variance  due  to  the 
noise  level  being  below  threshold.  As  the  %  jitter  is  increased  through  the  threshold  region  the 
estimates  first  show  occasionally  poor  results,  as  with  20%  jitter  in  figure  5.  These  statistically 
unusual  events  correspond  to  excessive  outlier  contamination  of  S'.  For  25%  jitter,  shown  in 
figure  6,  now  above  threshold,  many  more  values  of  T'  are  not  close  to  optimal.  The  situation 
increases  more  for  30%  jitter.  For  yet  higher  %  jitter,  values  of  T'  are  statistically  alike  with 
variance  comparable  to  TmeAi  as  shown  in  figure  4.  The  behavior  depicted  in  figures  5  and  6 
can  be  exploited  to  improve  T1  if  multiple  realizations  are  available,  as  we  consider  next. 

Multiple  Realizations 

If  multiple  independent  realizations  of  S  are  available  then  significant  improvement  in  T' 
can  be  obtained  in  the  transition  region  below  the  threshold.  Let  be  the  ith  estimate  based 
on  the  ith  realization  of  S.  Because  of  the  behavior  of  T‘  due  to  outliers  described  above,  we 
will  eliminate  values  of  that  are  excessively  contaminated  by  outliers,  and  then  average 
over  the  remaining  T^y  This  is  accomplished  by  rejecting  values  of  that  are  significantly 
distant  from  a  robust  measure  of  std{T^},  where  the  set  { X contains  the  estimates  of  T 
over  all  realizations. 
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Suppose  there  are  Nr  realizations  of  S.  We  form  a  robust  estimate  of  the  standard 


deviation  as  (see  Ljung  [16,  p.400]) 

*•  =  <46> 
where  a  (the  median  absolute  deviation)  is  the  median  of  {|T^  -  f  |}  and  T  is  the  median 
°f  Values  of  >  3 ar,  say,  are  rejected  as  outlier  contaminated,  leaving  N*r  <  Nr 

estimates  of  T .  We  then  average  over  the  remaining  N'r  values  of  to  yield  the  final  estimate 


1  N* 

:%y 


(47) 


Figures  7  and  8  demonstrate  the  improvement  possible  with  this  scheme.  Results  in  figure 
7  are  based  on  A  =  0.25,  x0  —  15,  and  15%  outliers.  Shown  are  normalized  MSE  of  Tmea 
and  T',  and  the  CRB,  all  based  on  a  single  realization.  Note  that  with  this  level  of  missing 
observations  and  contamination  T'  does  not  approach  the  CRB  until  the  %  jitter  is  less  than 
5%.  The  variability  in  the  MSE  of  Tmea  is  due  to  the  relatively  small  number  of  Monte  Carlo 
trials  being  inadequate  to  fully  characterize  all  possible  outlier  realizations. 

Also  shown  in  figure  7  is  Tm  for  Nr  =  5  realizations,  averaged  over  200  Monte  Carlo 
runs.  This  approach  yields  good  estimates  of  T  with  a  threshold  around  20%  jitter,  due  to  the 
judicious  elimination  of  when  it  is  strongly  outlier  contaminated.  Typically  2  or  3  of  the 
set  {T^}^-5  were  rejected.  These  results,  with  A  =  0.25  and  15%  outliers,  are  similar  to  those 
of  figure  1  (A  =  0.25,  no  outliers)  and  figure  4  (A  =  0.25,  5%  outliers).  Note  that  Tm ,  based 
on  multiple  records,  can  exceed  the  single  record  CRB  shown  in  figure  7.  Of  course,  the  price 
paid  in  using  Tm  to  combat  outliers  is  the  necessity  for  multiple  data  records.  This  is  to  be 
contrasted  with  simply  averaging  over  all  Nr  realizations  without  eliminating  any  of  the  T^ , 
which  exactly  corresponds  to  the  performance  curve  for  T'  in  figure  7,  obtained  by  averaging 
over  200  realizations.  Superior  performance  is  obtained  with  Tm  using  only  five  (versus  200) 
realizations. 

A  final  Monte  Carlo  experiment  is  shown  in  figure  8,  repeating  the  experiment  used  to 
generate  figure  7  but  now  using  A  =  0.50,  x0  =  10,  and  15%  outliers.  Here,  with  a  very  high  level 
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of  missing  observations  and  outliers,  T'  has  failed  to  achieve  the  CRB  for  even  small  percentages 
of  jitter.  However,  Tm  (again  with  Nr  =  5  realizations)  shows  good  performance  for  up  to  20% 
jitter.  In  this  difficult  case,  when  a  single  realization  is  insufficient  to  gain  improvement  over 
Tmea  with  the  refinement  algorithm,  multiple  realizations  can  be  successfully  exploited  to 
obtain  good  results.  It  can  also  be  noted  that  the  MEA  continues  to  yield  a  decent  estimate, 
based  on  a  single  realization.  Numerical  experiments  for  the  MEA  with  values  up  to  A  =  0.80 
and  10%  outliers  are  described  in  [6]. 

Remark 

1.  Estimation  of  (j)  and  T  is  a  linear  regression  problem.  For  example,  if  kj  =  N  +  1  —  j, 
j  —  1,2, ...  ,7V  then  the  problem  is  analogous  to  line  fitting  with  </>  playing  the  role  of 
intercept  and  T  the  role  of  slope.  The  observed  noisy  data  are  much  more  sensitive  to 
changes  in  the  slope  than  in  the  intercept,  making  estimation  of  the  intercept  considerably 
more  difficult.  The  CRB  of  the  ML  estimate  of  </>  is  0(AT-1),  while  the  CRB  for  ML 
estimation  of  T  is  0(N~3),  e.g.,  see  Kay  [13,  p.43].  The  estimate  <(>x  is  not  ML  and  has 
high  variance  for  even  moderate  noise  levels.  Thus,  if  it  is  desired  to  estimate  T  only,  it 
may  be  advisable  to  work  with  yj  =  tj  —  tj+i,  eliminating  the  need  to  estimate  (f). 


6  The  Common  Oscillator  Problem 


Many  radar  and  communications  systems  rely  on  a  very  stable  oscillator  to  provide  an 
accurate  time  reference.  Pulse  rates  will  then  generally  be  multiples  of  the  common  oscillator 
rate,  obtained  via  countdown  circuitry,  e.g.,  see  Wiley  [25,  sect.  8.6.2].  Examples  include  radar 
PRI  switching  and  pseudo-randomly  jittered  hop  times  in  frequency  hopping  communications. 
In  such  cases  it  may  be  desired  to  estimate  the  underlying  period  of  the  common  oscillator. 
This  can  be  accomplished  using  the  MEA  approach  provided  the  observations  are  of  the  form 
(1). 
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One  such  case  occurs  when  successive  pulse  intervals  are  multiples  of  the  common  oscil¬ 
lator  period.  Now,  (1)  is  applicable  with  %  —  fcj+i  =  £T,  where  £  is  a  positive  integer.  For 
example,  £  might  be  generated  pseudorandomly.  Thus  the  MEA  will  yield  an  estimate  of  T 
and  the  refinement  algorithm  may  be  used  to  enhance  the  MEA  estimate. 

7  Conclusions 

The  MEA  provides  an  efficient  tool  for  finding  the  period  T  from  the  data  set  S  arising 
from  (1).  It  is  especially  useful  with  many  missing  observations  and  contaminated  data.  With 
missing  observations  it  is  not  straightforward  to  apply  the  least-squares  minimum  variance 
solution,  and  the  addition  of  outliers  makes  the  problem  that  much  worse.  Without  prior 
knowledge  of  their  distribution  an  outlier  removal  strategy  is  not  obvious.  Because  the  MEA 
is  able  to  work  robustly  in  this  situation  it  may  also  be  exploited  to  provide  a  subset  S'  C  S 
that  may  be  used  to  enhance  the  MEA  period  estimate.  Although  not  in  general  outlier-free, 
S'  is  outlier  reduced. 

Nonlinear  estimates  of  the  kj' s  are  used  together  with  the  subset  S'  in  the  least-squares 
solution,  yielding  an  estimate  of  T  that  is  close  to  or  attains  the  CRB  for  moderate  noise 
levels  and  is  in  general  an  improvement  over  the  MEA  estimate.  The  enhanced  estimate  of 
T  exhibits  a  threshold  effect  due  to  the  nonlinearity  in  estimating  the  %’s,  with  estimates 
reaching  the  CRB  below  the  threshold.  In  the  case  of  contaminated  data  it  is  possible  to 
extend  the  threshold  to  higher  noise  levels  if  multiple  independent  data  records  are  available. 
This  cannot  be  accomplished  by  simply  averaging  over  records,  but  instead  relies  on  deleting 
those  estimates  of  T  that  are  highly  contaminated.  This  is  achieved  by  robustly  estimating  the 
standard  deviation  of  the  estimates  over  records,  and  eliminating  those  that  are  of  significantly 
high  statistical  variation. 

Further  improvement  is  likely  to  be  achieved  if  sort  parameters  other  than  simply  the  time 
of  arrival  can  be  incorporated,  which  is  a  topic  for  further  research.  When  dealing  with  more 
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complex  problems  such  as  radar  PRI  estimation  in  the  presence  of  random  period  changes, 
it  is  likely  that  parameters  such  as  direction  of  arrival,  pulse  duration,  pulse  amplitude,  etc., 
will  be  necessary,  e.g.,  [17,18].  Also  interesting  might  be  the  combination  of  the  MEA  with 
deinterleaving  problems.  For  example,  a  period  may  be  postulated  by  time  differencing  arrival 
times,  and  then  selected  data  points  within  a  window  analyzed.  In  this  sense  the  MEA  provides 
another  tool  to  use  in  combination  with  existing  methods. 


8  Appendix 


Assume  a  noise  pdf  of  the  form 


-a  cos(27T77/T)  rp 

=  2 */„(*)  '  W  S  2  • 


(48) 


Here,  /(r /)  =  f(—rj)  for  all  a  implies  zero-mean.  Assuming  that  (3  =  [<j>,  T]T  is  known,  the  pdf 
of  tj  given  (3  is 

f&\0)  =  m  -  ksT  -  4>),  (49) 

and  for  iid  noise, 


N 


Now  the  log-likelihood  function  is 


N 


2? r, 


L((3)  =  ln/(t|/3)  =  X)[acos^r(^  “  -ln2irl0(a)]. 


3= 1 


So,  the  ML  estimate  is  found  by  maximizing  L((3)  with  respect  to  F,  yielding 


^  /  (ti-tf))' 

TMl  =  argmaxr  cos  ^2?r - — - 


Using  the  fact  that 


cos 


then,  for  (f>  constant,  maximizing  L((3)  is  equivalent  to 

Tx  =  argmaxy 

which  is  in  the  form  of  spectrum  analysis  of  a  point  process  with  occurrence  times  tj. 


N  t- 

3= 1 


(50) 


(51) 


(52) 


(53) 


(54) 
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Table  1:  Comparison  of  data  from  figure  1. 


%  jitter 

Tmea  mean(std) 

- 5*; - - - 

T  mean(std) 

VGRB 

2 

1.0000  (0.0007) 

1.0000  (1.826  x  10"6) 

1.829  x  10"5 

4 

1.0000  (0.0011) 

1.0000  (3.587  x  10"5) 

3.669  x  10~5 

6 

1.0003  (0.0016) 

1.0000  (5.056  x  10“5) 

5.532  x  10"5 

8 

1.0000  (0.0022) 

1.0000  (7.427  x  10~5) 

7.366  x  10"5 

10 

1.0002  (0.0028) 

1.0000  (9.275  x  10"5) 

9.228  x  10~5 

15 

0.9998  (0.0043) 

1.0000  (1.438  x  10“4) 

1.387  x  10-4 

20 

1.0000  (0.0054) 

1.0000  (1.907  x  10"4) 

1.836  x  10“4 

25 

0.9999  (0.0072) 

1.0000  (2.210  x  10"4) 

2.308  x  10-4 

30 

1.0007  (0.0084) 

1.0000  (1.583  x  10"3) 

2.758  x  10~4 

35 

0.9999  (0.0110) 

0.9994  (3.401  x  10~3) 

3.222  x  10"4 

40 

1.0008  (0.0139) 

0.9996  (6.216  x  10"3) 

3.672  x  10“4 

45 

1.0052  (0.0426) 

0.9999  (3.285  x  10-2) 

4.178  x  10-4 

50 

1.0133  (0.0741) 

1.0041  (5.830  x  10-2) 

4.792  x  10-4 

Figure  1:  Monte  Carlo  estimation  results  (without  outliers,  A  =  0.25,  N  =  100,  xq  =  5).  Solid 
=  CRB,  x  =  TmeAi  dash  (*)  =  T,  and  dash-dot  (-f )  =  T  incorporating  <j)x. 
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Figure  2:  Monte  Carlo  estimation  results  (without  outliers,  A  =  0.50,  N  =  100,  cc0  =  5).  Solid 
--  CRB,  x  =  Tmea ,  and  dash  (*)  =  T. 


Figure  3:  Monte  Carlo  estimation  results  (without  outliers,  A  =  0.75,  N  —  100,  x$  —  5).  Solid 
=  CRB,  x  =  Tmea ,  and  dash  (*)  =  T. 


Figure  4:  Monte  Carlo  estimation  results  (with  5%  outliers,  A  =  0.25,  N  =  100,  x0  =  15). 
Comparison  of  T  versus  T'  (full  versus  selected  data  reuse  in  the  refinement  algorithm).  Solid 
=  CRB,  x  =  Tmea j  dash-dot  (o)  =  T,  dash  (*)  =  T'. 


Figure  5:  Realizations  of  T'  for  different  percentages  of  jitter  near  threshold,  with  experimental 
parameters  matching  those  of  figure  4. 
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Figure  6:  Realizations  of  T'  for  different  percentages  of  jitter  above  threshold,  with  experimen¬ 
tal  parameters  matching  those  of  figure  4. 


Figure  7:  Multiple  realization  results.  Solid  =  CRB,  x  =  Tmea ,  and  dash  (*)  =  T',  for  single 
data  records  (A  =  0.25,  N  =  100,  xq  =  10,  and  15%  outliers).  Also,  dash-dot  (o)  =  Tm  (Nr  —  5 
data  records). 
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Figure  8:  Multiple  realization  results.  Solid  =  CRB,  x  =  TmeAi  and  dash  (*)  =  T',  for  single 
data  records  (A  =  0.5,  N  —  100,  xq  =  10,  and  15%  outliers).  Also,  dash-dot  (o)  =  Tm  (Nr  =  5 
data  records). 
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